Orbit Determination and Prediction, 
and Computer Programs 

By A. J. CLAUS, R. B. BLACKMAN, E. G. HALLINE 
and W. C. RIDGWAY, III 

(Manuscript received March 15, 1963) 

Orbit determination and prediction programs are needed to generate 
ephemerides for the satellite. Orbit determination is from tracking data 
consisting of angles only, and is based on a modified version of a method 
by R. E. Briggs and J. W. Slowey of the Smithsonian Institution. Trends in 
the data due to perturbations from a Keplerian orbit are removed before 
this process, and estimates of the orbital elements from individual passes 
are combined statistically to produce refined estimates. Ephemeris calcu- 
lation is by a semi-analytic method in which deviations from a Keplerian 
orbit are obtained by integrating the perturbing forces. The programs to 
implement these procedures have been written for both the IBM 7090 and 
the IBM 1620 computers. 

I. INTRODUCTION 

The following paper describes the methods and programs used in the 
Tclstar project for the purposes of orbit determination and ephemeris 
calculation. The orbit determination process involves the computation 
of orbital elements from tracking data obtained during each pass, and 
subsequent refinement by combining such single-pass estimates. The 
tracking data arc in terms of angular observations only. The ephemeris 
calculations involve standard procedures for computation of Keplerian 
orbits and perturbations due to the earth's oblateness. 

It is well known that in the problem of orbit determination from 
angular data only, three observations (each observation consisting of two 
angles and a time) are not sufficient to determine an orbit if the three 
sightlincs are. coplanar. If the three sightlines are nearly coplanar, the 
computed orbital elements may reflect large uncertainties which are not 
necessarily due to observational errors. Hence, the method used is based 
on the determination of a set of orbital elements from four observations. 
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This method is a modified form of the method described in Ref. 1. In 
the modified form, initial estimates of the orbital elements are computed 
from the first and last of the four observations, supplemented by esti- 
mates of the ranges corresponding to them (Section 2.1). By a method 
of successive approximations the two ranges are adjusted until an agree- 
ment with the second and third observations is secured in a least-squares 
sense (Section 2.2). 

With typically many more than four observations in one pass through 
the visibility zone of an angular tracker, the observations are divided 
into four nonoverlapping blocks, each block containing the same number 
of observations, N. Taking one observation at a time from each block, in 
serial order, N sets of orbital elements are computed. These N sets are 
combined into a single set of intrapass average orbital elements and an 
associated covariance matrix (Section III). 

Trends in the data due to perturbations induced by the earth's ob- 
lateness are removed by a method which is essentially the same as that 
described in Ref. 2 (Section IV). 

Sets of intrapass average orbital elements, and their associated 
covariance matrices, from two or more passes, are combined into a set 
of interpass average orbital elements and an associated covariance 
matrix (Section V) . The method used is similar to the method described 
in Refs. 3 and 4. inasmuch as it was motivated by the desire to avoid 
the necessity of pooling all of the observational data from two or more 
passes in order to derive refined estimates of the orbital elements, as 
would have to be done in the classical "differential corrections" method 
commonly ascribed to K. F. Gauss (1777-1855). However, the method 
used differs from the referenced method in two respects, viz., (a) the 
covariance matrix associated with each set of intrapass average orbital 
elements is related to the actual observational data for the pass, and 
(b) the necessity of computing the partial derivatives of all of the ob- 
served angles (numbering 8N in each pass) with respect to each of the 
orbital elements is avoided. On the other hand, this method gives single- 
pass estimates of the orbital elements which are biased even when the 
observational errors are not biased. These biases may be appreciable for 
short passes associated with low altitudes of the satellite near perigee. 
Methods for removing or reducing these biases have been under study 
but were not ready for use before the launching of the Telstar satellite 
on July 10, 1962. 

This orbit determination method was designed to permit effective 
antenna pointing operations with the use of a modest computing facility. 
The program implementation (Section VI) consists of two major 
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program subsystems. The first of these is the orbit determination pro- 
gram (Section VIII), which determines the characteristics of the satellite 
orbit from tracking information. The second program (Section IX) com- 
putes orbit predictions from a knowledge of these orbit characteristics. 
The operational results obtained in using these methods and programs 
are discussed (Section X). 

II. ORBIT DETERMINATION FROM ANGLE-ONLY DATA 

2.1 Orbit Determination from Two Observations and Estimated Ranges 

Two observations (each observation consisting of two angles and a 
time) and estimates of the ranges (i.e., topocentric distances) along the 
two sightlines are sufficient to establish two points P x and P 2 through 
which an orbit can be passed at the times of observation by only one 
set of orbital elements. Denoting the geocentric distances by /•] and r 2 , 
and the geocentric angular difference by i2 , we have 



_ ail — I 2 — m~) 

ri r+r~ ~ • 

a{\ — I — m 



(0 
(2) 



(3) 



1 + ' cos 0i2 + m sin d vl ' 

where 

/ = e cos co, 

m = e sin to, 

a is the semi-major axis, c is the eccentricity, and to is the argument of 
perigee referred to the first sightline. From (1) and (2), we have 

ail + a-itn = a? (4) 

where 

a n 
ay = COS Pi2 , 

r-i 
a 2 = sin 0i2 , 
?"i — r-i 



a 3 = 



r-> 



It is convenient to regard either I or m as an independent variable. 
Actually, in order to avoid an indeterminacy and to improve accuracy, 
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preference is given to Z if | «i | ^ I «2 |, and to m if J a% \ > | a 2 1. In either 
case, the other two of the three quantities I, m, and a are determined by 
(1) and (4). These, in turn, through Kepler's equation, determine a 
travel time between Pi and P 2 • A search is then made for the value 
of I which gives the observed travel time. For bounded orbits (the 
only ones of interest for ground-to-ground communications) the search 
is confined to the interval (/_ , l+), where 

J± = T~i * («i«3 ± | "2 1 Vai 2 + a? — tt3 2 )- 

ai 2 + at 

With a, I, m (hence also e and u) determined, the time of perigee pas- 
sage, t, is determined through Kepler's equation. 

2.2 Orbit Determination from Four Observations 

The two observations involved in the procedure described in the 
preceding section are the first and last of a set of four observations. 
Subsequent to that procedure, the four angles corresponding to the 
times of the second and third observations are computed and compared 
with the observed angles. The sum, $, of the squares of the differences 
between the computed and observed angles is regarded as a function of 
the estimated ranges D\ and Z) 4 associated with the first and last of the 
four observations. The quantity <£ is next minimized with respect to 
Di and Z) 4 by a method which is analogous to the classical "differential 
corrections" method. (With only three observations corresponding to 
coplanar sightlines there would be only one angular difference, and 
therefore Di and Da would be indeterminate.) This method involves the 
solution of two simultaneous equations which are linear in the correc- 
tions to D\ and D 4 , with coefficients which are quadratic in the first- 
order partial derivatives of the computed angles with respect to D Y and 
Da . The terms which do not involve the corrections to Di and 7) 4 arc 
products of the first-order partial derivatives and the angular differ- 
ences. Since the partial derivatives are functions of Di and D 4 , the 
minimization of $ is an iterative procedure which is terminated when 
the values of a, c, w, and r are sufficiently stabilized. Detailed formulas 
are given in Ref. 5. 

With a, e, w, and r determined, the orientation of the plane through 
Pi , P 2 , and the center of the earth gives the values of £2 and i, where 
fi is the longitude of the ascending node and i is the inclination of the 
orbital plane. 
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III. ORBIT DETERMINATION FROM 4N OBSERVATIONS IN ONE PASS 

The combined procedures described in Sections 2.1 and 2.2 are applied 
to as many sets of four observations as may be drawn from all of the 
reliable observations for each pass in accordance with the method of 
selection indicated in the third paragraph of Section I. The N sets of 
orbital elements are then combined into a single set of intrapass average 
orbital elements. In addition, an associated 'co variance matrix (an es- 
timate of the variability of the mean in a sample of size iV drawn from 
a correlated multivariate population) is computed in accordance with 
the standard formulas 



^2 («»• — «)" 



fc = ££ 



N{N - 1) ' 



E («< - a) (ft - 0) 

C ^ = ' N(N - 1) 



where a stands for each of the six orbital elements (with average a), 
and for each of the other five (with average 0). 

A typical result of the single-pass routine, as described up to this 
point, is shown in Tables I and II. The orbital elements listed as "exact 
value" were used to generate tracking angles. These angles, combined 
with random errors from a normal population with a standard deviation 
of 0.2 milliradian, were processed. It may be noted that were it not for 
the strong correlation between some of the orbital elements, errors in 

Table I 



Exact value 
Sample mean 
Standard devia- 
tion of sample 
mean 



S2, degrees «", degrees a, feet 



144.4462 

144.4455 

0.0018 



46.9190 

46.9184 
0.0018 



31,507,194 0.240764 

31,573,342 0.240688 

6,212 0.000120 



0), degrees r, seconds 



171.6756 

171.6349 

0.0374 



47,953.227 

47,950.120 

2.852 



Table II — Correlation Coefficients 
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1 

0.13189 

0.64082 

-0.57898 

-0.67162 

-0.64998 


0.13189 
1 
-0.40292 

0.52030 
0.44004 
0.47200 


0.04082 
-0.49292 

1 
-0.90559 
-0.99058 
-0.99680 


-0.57898 
0.52030 

-0.96559 
1 
0.92798 
0.94376 


-0.67162 
0.44004 

-0.99058 
0.92798 

1 
0.99797 


-0.64998 
0.47200 

-0.99680 
0.94370 
0.99797 
1 
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the elements of the order of the standard deviations shown could result 
in pointing angles for the same pass with errors as large as 10 times the 
standard deviation of the original tracking errors. 

A caveat should also be noted with respect to the precision of compu- 
tation of the covariance matrix. Any matrix which purports to be a 
covariance matrix must have a nonnegative determinant. Due to the 
high correlation among the elements a, e, w, and r, however, values of 
10~ for the determinant of the correlation matrix are common. Errors 
of the order of 0.1 per cent in some covariances could result in a matrix 
with a negative determinant. Such a matrix can still serve as a guide in 
judging the reliability of the orbital elements obtained, but the use of 
this matrix for interpass orbit refinement would very likely lead to 
absurd results, such as negative variances. 

IV. TREND REMOVAL 

Since the procedures described in Sections 2.1 and 2.2 are based on 
the assumption that the orbit is Keplerian, it is important to determine 
the extent to which it is necessary and sufficient to correct for deviations 
from that assumption. Such deviations, usually called perturbations, are 
induced by the asphericity of the earth, drag, radiation pressure, etc. 
Preliminary computations, confirmed by tests with artificial data, in- 
dicated that for the orbit and satellite under consideration here it 
would be necessary and sufficient to correct only for the earth's 
oblateness. The corrections are made to the observational data. Detailed 
formulas for the corrections are given in Ref. 5. These formulas involve 
the orbital elements which, however, do not need to be known to high 
accuracy for the purposes of trend removal. If sufficiently accurate 
values of the orbital elements are not available for trend removal, they 
may be obtained by including trend removal in the iterative routine of 
Sections 2.1 and 2.2 after the first values of the orbital elements have 
been obtained without trend removal. 

Table III shows the importance of trend removal for the effects of 
oblateness. The same input data, which included the effects of oblate- 
ness, were used in both runs. The errors in the second run (without trend 
removal) are not acceptable. In particular, the error in the semi-major 
axis could lead to an error in predicted pointing angles of as much as 
1.5° after only one period. 

Table IV shows the speed of convergence, with trend removal, in the 
absence of initial estimates of the orbital elements. After only one 
iteration (one-half minute additional computing time for 200 observa- 
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Table III 



Exact osculating 
elements at the 

center of the pass 
Results of run no. 1 
(with trend re- 
moval) 
Results of run no. 2 
(without trend re- 
moval) 



S!, degrees 



144.4439 
144.4439 



i, degrees 

4G.9170 
46.9169 



a, feet 



144.4372 40.9152 



31,566,742 
31,560,884 
31,542,821 



0.240879 
0.240875 
0.241313 



o, degrees 



171.6124 
171.6118 

171.75(H) 



t, seconds 



47,950.421 
47,950.365 
47,961.173 



Table IV 



Exact osculating 
elements at the 
center of the pass 
Results of run no. 1 
Results of run no. 2 
Results of run no. 3 



ft, degrees 



144.4439 



144.4377 
144.4439 
144.4439 



i, degrees 



46.9170 



46.9149 
46.9169 
46.9169 



a, feet 



31,566,742 



31,543,423 
31,566,909 
31,566,886 



0.240879 



0.241311 
0.240874 
0.240875 



w, degrees 



171.6124 



171.7487 
171.6117 
171.6118 



t, seconds 



47,950.121 



47,960.978 
47,950.352 
47,950.302 



tions, on an IBM-7090 computer), acceptable orbital elements were 
obtained. 

V. COMBINATION OF SINGLE-PASS ORBITAL ELEMENTS 

The method of combining single-pass estimates of the orbital elements 
is based on a matrix formula derived briefly as follows. Let .? be a vector 
(.i.e. a one-column matrix) estimate of the vector z, with ave [x — z\ =0 
and cov {.? - z\ = A, where .4 is a covariance matrix. Similarly, let 
y be another estimate of z, with ave [y - z\ = and cov \y — z] = B. 
If x and y obey independent multivariate normal probability distribu- 
tions, the "maximum likelihood" estimate of z is the z which minimizes 
the quadratic form 



Q = (z - ty-A- l -(i - ■(•) 
where the primes denote transposition. Thus, 



z - yY-B 



(z - y) 



A l -(z - x) + B~ l -{z - y) = 0. 



whence, 



z = (A 



1 + B' 1 )- 1 - 



(A~ l x + B~ l y) 



(5) 
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with covariance matrix 

C = (A~ l + B' 1 )- 1 - (6) 

In fact, it may be easily verified that 

Q = [z - C-(A~ l x + BT l y)\'-Cr l '[2 - C-{A~ x x + B"*f)] 

+ terms independent of z. 

A somewhat longer derivation without the normality assumption, in 
which the main diagonal (variance) elements of C are minimized, leads 
to the same results. 

Formulas (5) and (G) require three matrix inversions which result 
in an intolerable loss of accuracy in cases of highly correlated estimates 
of the orbital elements. This difficulty is relieved to a very large extent 
by using the equivalent formulas 

z = wix + w 2 y - (wiP - w*Q)(P + Q)~\x - y), (7) 

C = \[w x A + wo_B - ( Wl P - w 2 Q)(P + Q)~\A - B)] (8) 

where P = AG, Q = BG, G is an arbitrary six-by-six matrix, and w\ , 
w 2 are any two six by six matrices whose sum is a unity matrix (see 
Appendix A). Formulas (7) and (8) require only one matrix inversion. 
The matrix G can be constructed so that the matrix (P + Q) is well 
suited for inversion. 

As a matter of additional necessity, formulas (7) and (8) were 
further transformed by the introduction of matrices U, V, denned by 
U = SA S, V = SBS, where S is a diagonal matrix whose elements are 

Sa = (An + Bn) , 

so that the diagonal elements of the matrix ( U + V) are unity. Re- 
stricting Wi , w 2 to diagonal matrices, then, 

z = Wix + w 2 y — R(x - y), (9) 

C = \[wxA + w 2 B - R(A - B)], (10) 

where 

R = S-\ Wl P - w 2 Q)(P + Qy'S, (11) 

P = UH, Q = VH, and H = S~ l G. The formal construction of the 
arbitrary matrix H is not necessary. The matrices P, Q, and P + Q 
are obtained by linear combinations of rows and/or of columns of the 
matrices U, V, and U + V according to rules which are easily pro- 
grammed for a digital computer. 

Two details must be noted in the use of these formulas for combining 
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sets of orbital elements. The first detail is that the orbital elements 
are actually "oscillatory" orbital elements which vary with time; 
therefore, each set is necessarily referred to a specific "epoch." Hence, 
before combining two sets, the set referred to the earlier epoch must 
be "updated" to the later epoch. In updating a set of orbital elements, 
it must also be noted that the "time of perigee passage" is actually 
the "time of with perigee passage," where m is a specific number, 
usually different from the one for the set referred to the later epoch. 
The second detail to be noted is that the covariance matrix for the set 
referred to the earlier epoch must also be updated. 

If C'i is the covariance matrix to be updated, the updated covariance 
matrix is given by the formula 

C'i = JCiJ', 

where J is the Jacobian of the updated orbital elements with respect 
to the orbital elements from which they were predicted. Even in the 
hypothetical case of Keplerian orbits, in which all of the orbital elements, 
with the possible exception of t, arc constants, the Jacobian may differ 
from a unity matrix. For example, if the updating is through m times 
the period 2tt\/ a^/k, so that 

t\ = fi + 2irm\/a i 3 /l:, 

then, 

dn/dai = ZTmy/di/k. 

The results of a test problem of this hypothetical sort are shown in 
Table V, in which the updating was through one period. The standard 
deviation of the improved estimate of the semi-major axis is approx- 
imately i4o times the average of the corresponding standard deviations 
for the two runs. The improved estimate is in error by only 52 feet. 
Table VI shows the results of a more realistic test problem in which 
the input data included perturbations due to the earth's oblateness. 
With "no updating" of the orbital elements and the covariance matrix 
from the earlier pass, except only to the extent required in the hj r po- 
thetical case of Keplerian orbits, the "improved" semi-major axis is in 
error by 5094 feet, which is inconsistent with the standard deviation of 
only 73 feet. However, with updating of the orbital elements, taking 
account of the effects of the earth's oblateness, the error is only 72 feet, 

VI. PROGRAM DESCRIPTION 

The computer program system required to track a satellite and gener- 
ate steering information for the communications antenna is divided into 
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two major subsystems. The first of these is the orbit determination 
program, TELETRACK, which determines the characteristics of the 
satellite orbit from tracking information. The second major program, 
TELEPATH, computes orbit predictions from a knowledge of these 
orbital characteristics. 

The division of the program system into these two parts is not only 
natural, but is also dictated by systems considerations. One of the 
requirements on the system was to minimize the amount of data trans- 
missions. Ephemeris data to steer the communications antenna can be 
generated from the six orbital elements, and a division of the program 
system into two components linked together only by these six numbers 
achieves this requirement if each ground station is provided with suitable 
computational facilities. Stations having communications antennas 
require the program TELEPATH and updated sets of the orbital 
parameters. Stations having tracking antennas process the tracking data 
with TELETRACK and broadcast the updated elements to other 
stations as they become available. 

The IBM 1620 computer was chosen to provide on-site computations. 
The IBM 7090 computer was used, however, for the initial development 
of the program systems. This was done for two reasons. First of all it 
was desirable to take advantage of the more powerful facilities and speed 
of the larger computer to facilitate the development and testing of the 
methods employed in the program system. Secondly, it was desirable 
to have the complete program system available at the Whippany, N. J., 
location of Bell Telephone Laboratories as a back-up to the on-site 
computer centers. Experience has shown that it is absolutely essential 
to have these duplicate programs available for testing and checking of 
the on-site operations. 

By the nature of the 7090 and 1620 computers, different operating 
philosophies are required for each. The speed of the 7090 and turn- 
around times inherent in a large computation center are such that the 
programs must be as automatic as possible. However, they must also be 
flexible enough to allow selected programs from the system to be per- 
formed when necessary. Towards this end the following system evolved. 
The entire set of 7090 programs can be run consecutively as a single 
automatic chain job. Each program communicates to the following 
program through a magnetic tape, but as far as the computation center 
is concerned each program is a separate job. As a consequence, each 
program can also be run independently (with input provided by cards) 
since it is an entity in itself. The hidden gain in this system is the fact 
that there is only the one flexible version of each program, thus elimi- 
nating confusion and mistakes. 
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For the 1620, which is devoted entirely to this problem, and which is 
a slower machine, such completely automatic operation is not necessary. 
The system can be run automatically, but is usually run with more 
direct operator intervention. This allows greater flexibility and the 
ability to monitor intermediate results. On the 1620 the two major 
program systems are broken down conveniently into several program 
components. Each of these programs runs independently of the others, 
receiving input data generated by one of them and preparing output 
data for another. Operation of the program systems is achieved by 
loading and running one of the program components at a time. The 
various program components are stored on magnetic tape, and each 
program in the system loads the next program into the computer from 
this tape. Transfer of data between the programs is accomplished by 
punched cards, magnetic tapes and common memory storage. The 
method of data transfer in a particular instance depends upon the nature 
and quantity of the data. 

Numerous error conditions were anticipated while the programs were 
being written. Many of these are handled automatically by the programs 
themselves. Some must be taken care of by manual intervention. 

VII. INERTIAL COORDINATES AND ORBITAL ELEMENTS 

All orbital calculations must, of course, be referenced to an inertial 
(or near inertial) coordinate system. The basic system used in these 
programs is the usual earth-centered, right-handed rectangular system. 
The X-Y plane coincides with the earth's equatorial plane, the X-axis 
is parallel with the line of equinoxes, and the Z-axis passes through the 
North pole. The orientation of the earth in this system at the time of an 
observation is obtained from UT 2 at time of observation and the Green- 
wich Mean Sidereal Time at hours UT of date. Conversion from Mean 
Sidereal Time to Apparent Sidereal Time is made using the Equation 
of Equinoxes at hours UT of date; interpolation of this number to the 
time of observation was deemed unnecessary. 

The satellite orbit is described by means of the osculating orbital 
elements, consisting of 

(a) semi-major axis 

(b) eccentricity 

(c) right ascension of ascending mode 

(d) inclination angle 

(e) argument of perigee, and 

(f) time of perigee passage. 
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These elements specify the ellipse oscillatory to the satellite orbit at 
some instant in time. These six numbers are therefore accompanied by 
an epoch specifying the time of osculation. The time of perigee passage 
specifies the perigee passage immediately preceding the epoch and is 
stated in seconds relative to the epoch. 

The following paragraphs describe in some detail the two program 
systems, TELETRACK and TELEPATH. 

VIII. TELETRACK PROGRAM SYSTEM 

The TELETRACK program system processes tracking data in terms 
of azimuth and elevation to produce estimates of the six orbital elements 
describing the satellite orbit. It processes tracking data from one pass 
over the tracking station at a time to produce a "single-pass estimate." 
Single-pass estimates are combined to provide "combined estimates." 
The combining of several single-pass estimates provides a statistical 
averaging of the several independent estimates and a refinement based 
on the separation in time of the various independent estimates. 

A flow chart of TELETRACK is shown in Fig. 1. Each of the major 
program components and the modes of data transfer between them arc 
shown. A few of the program switches which control the mode of opera- 
tion of the system are also shown. 

8.1 TELED 

TELED is the input/edit section of TELETRACK. Inputs to this 

program are 

(a) tracking data consisting of time, azimuth and elevation for one 

pass, and 

(b) data cards containing date and number of pass, identification of 
the tracking station and satellite, meteorological conditions during the 
pass, GMST at (V 1 of date, estimates of the orbital elements, number of 
data sets to be selected {N), and values of the mode control switches 
for TELETRACK. 

TELED reads the tracking data from tape and performs format and 
units conversion. Data points for which the precision tracker was not in 
autotrack or for which the signal-to-noise ratio level was not above a 
predetermined level (usually 4 or 5 db) are rejected. Furthermore, data 
points for which the elevation is below 7.5° or above 82.5° are rejected. 
The specified number (4N) of data points is selected from the group 
satisfying these criteria. The set of data so selected is distributed as 
uniformly as possible over the available set. 
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Fig. 1 — TELETRACK flow chart. 



Boresight and refraction corrections are then applied to the selected 
data points. Following this, coordinate conversions are performed to 
transform the data from the topocentric azimuth-elevation system to 
the inertial coordinate system. Deviations of the vertical from the 
normal to the geodetic spheroid are accounted for in this process. Since 
range data are not available, the results of the coordinate conversion are 
in terms of the direction cosines relative to the inertial system of the 
observed sight lines. Also computed are the coordinates in inertial space 
of the tracker at the time of each observation. 

Outputs from TELED are stored on magnetic tape for subsequent 
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programs. The primary output consists of 4iV data points. Each data 
point contains time of observation, the three direction cosines of the 
observed sight line, sine and cosine of the right ascension of the tracker 
at the time of the observation, and the inertial coordinates of the 
tracker at the time of the observation. Various other data are also 
stored on this tape for subsequent programs. 

8.2 TREND 

The direction cosines produced by TELED are adjusted by TREND 
to produce the set which would have been obtained had the .satellite 
been moving in an unperturbed, elliptical orbit throughout the pass. 
These adjustments are described in detail in Section IV above. The time 
of osculation (l c ) between the perturbed and unperturbed orbits was 
selected by TELED to correspond to the center of the pass and is 
passed on to TREND via tape. 

As noted above, an estimate of the orbital elements at time t B is 
needed. These are obtained by updating to time t c the elements sup- 
plied on the input cards to TELED. Program ORBFIX is used for this 
purpose, details of which are given below. 

The output from TREND consists primarily of the adjusted direction 
cosines for each observation. These are stored on a tape which is identical 
in format with the TELED output tape. By making these formats 
identical it is possible under one of the modes of operation to bypass 
TREND if estimates of the orbital parameters are not available. 

8.3 ORBFIX 

As mentioned above, ORBFIX updates a set of osculating orbital 
elements valid at one epoch to another epoch. The program essentially 
makes use of the subroutine OBLATE with only minor additional 
bookkeeping operations. The subroutine OBLATE is a numerical 
integration routine in true anomaly which integrates in steps of 0.08 
radian the first-order oblateness perturbation equations to provide the 
desired corrections. The equations also include sufficient second-order 
terms to allow taking steps of 2tt, so that in actual use steps of 2v are 
taken until a value w or closer to the desired point is reached. The 
program then integrates either forward or backward in small steps to 
reach the desired point exactly. It is also possible to go only in 2ir steps 
in cases where only limited accuracy is required. This results in a large 
time saving. 
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8.4 ORB EL 

Calculation of the orbital parameters is performed by the program 
ORBEL. The input data normally consist of the adjusted direction 
cosines from TREND. In the absence of initial estimates of the orbital 
elements, however, ORBEL can process the unadjusted data from 
TELED. The 4iV observations are divided into four nonoverlapping 
groups. A set of four observations is obtained by selecting one observa- 
tion from each group. N independent estimates of the orbital elements 
are calculated from the resultant N sets of observations. Averages, 
variances and covariances of the six elements for one pass are calculated 
from these. Details on the methods are given in Sections II and III. 

The estimates of the ranges required in producing the first set of 
elements are normally produced by TREND during the trend removal 
procedure. In the absence of trend removal these estimates must be 
supplied on the input data cards to TELED. Subsequent estimates of 
range are derived by ORBEL itself from its previous estimates of the 
elements. 

The output from ORBEL consists of a set of cards (an "ORBEL 
deck") containing the single-pass estimates of the orbital elements, 
the standard deviations of those elements, and the correlation coeffi- 
cient matrix. Pass number and the corresponding epoch are also stored 
on these cards. These cards are filed away for possible future use. 

The information on these cards is also retained in memory for use 
by the combination of passes program, COMPS. 

8.5 COMPS 

Combination of the estimates from the various passes is accomplished 
by the program COMPS. The method employed is described in Section 
V above. The inputs consist of two sets of orbital elements, standard 
deviations and correlation coefficient matrices. The first set, obtained 
from input cards, is either from a single ORBEL run or from an earlier 
COMPS run. The second set is from the current ORBEL run and is 
usually supplied directly by ORBEL through common memory storage. 
Under some modes of operation, however, the second set is supplied by 
cards. 

The output from COMPS is a set of cards (the "COMPS deck") 
identical in content and format with the ORBEL deck. These cards 
are filed to maintain a permanent record of the combined orbital ele- 
ments. The output data also replace the data from the first input set 
in memory in case certain operating modes are selected. 
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8.G Modes of Operation 



Several mode-control switches are provided to permit selection of 
one of a number of possible operating modes. The more significant of 
these switches are shown in Fig. 1. Each switch is identified by a num- 
ber and consists of a one-bit variable which is read from an input data 
card and stored in memory. 

In the normal mode of operation it is assumed that some estimate of 
the orbital elements is available for trend removal and that a combina- 
tion of the ORBEL output with an earlier COMPS output is wanted. 
All switches are set in the "off" condition and the sequence of opera- 
tions is TELED, TREND, ORBEL and COMPS in that order. The 
first set of inputs to COMPS is determined by the operator, who se- 
lects the proper COMPS deck, and the second set is supplied directly 
by ORBEL. 

An alternative mode of operation is to stop the program after the 
single-pass estimate is produced by ORBEL and then combine a number 
of such estimates in a "batch combination" at a later time. This is 
accomplished by turning switches 3 and 4 on and accumulating a 
number of ORBEL decks. These decks are fed to COMPS in order by 
time, with the earliest deck first. COMPS reads the first two decks 
and combines them, producing a combined estimate valid at the time 
of the second set. This in turn is combined with the third set to produce 
a combination of the first three decks valid at the time of the third. 
This process continues until all decks have been combined into a single 
estimate valid at the time of the last set. 

Another mode of operation is available in case estimates of the 
orbital elements are poor or unavailable. By turning switch 1 on, trend 
removal is skipped initially, and ORBEL is given unadjusted data 
with which to estimate the elements. If switch 2 is also on, ORBEL 
will call on TREND after computing this initial estimate of the orbital 
elements. This estimate is passed on to TREND for use in adjusting 
the data. Switch 2 is turned off, the data are adjusted, and then ORBEL 
is called upon a second time, this time to process data with trend clue 
to perturbations removed. 

IX. TKLEPATH PROGRAM SYSTEM 

The ephemeris generation for the Telstar satellite is carried out by a 
trio of programs collectively known as the TELEPATH program. 
The three individual programs are called MUVIS, COKE, and ACEXP, 
and are complete entities in themselves, solving distinctly separate 
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portions of the problem. The MUVIS program is solely concerned with 
finding times of future visibility or mutual visibility and updating the 
orbital elements to these time periods. Its output is a listing of future 
passes which is in itself useful, and a set of cards which serves as input 
to the COKE program. The COKE program generates a theoretical 
ephemeris for each pass as determined by the input cards, and outputs 
it on tape. The COKE program can also be used by itself to re-create 
any pass for which orbital elements are available. Both programs 
exist in almost identical form both for the IBM 7090 and the IBM 
1620. The only differences in the programs are due to storage limitations 
in the 1620. This results in some extra tape manipulations in the 1620 
programs which are unnecessary on the 7090. The final program, 
ACEXP, exists only on the 1620 and is used for adding predistortion 
and refraction corrections to the theoretical ephemeris. 

Fig. 2 shows the flow chart of the 1620 program with its various 
operating options. A more detailed description of the program follows, 
without reference to machine. 

9.1 MUVIS 

This program takes a set of osculating orbital elements at an epoch 
and using them predicts when the satellite will be visible at a designated 
site, and when it will be mutually visible with a second designated site. 
The emphasis in this program is speed with only a limited amount 
of accuracy. It is envisaged that this program will be used for planning 
and general information, and thus the methods used were chosen with 
this in mind. 

Basically, the program steps time by some increment, predicts the 
satellite's position in inertial coordinates for the new time, checks for 
visibility and mutual visibility, and continues. There is naturally a fair 
amount of bookkeeping associated with executing these steps, but they 
are essentially the heart of the program. 

Since the program consists of many iterations through the basic 
loop outlined above, it was felt worthwhile to streamline it as much as 
was possible. Towards this end the following steps were taken. 

(i) The program takes variable time steps. A coarse step is used 
until visibility is determined, and at this point a finer step is used for a 
more refined estimate. This feature is carried one step further by per- 
mitting a time step of close to a full period after visibility ends ; or when 
the satellite appears to be moving away from visibility. 

iii) When the satellite's position is calculated at some time, osculat- 
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Fig. 2 — TELEPATII flow chart. 

ing orbital elements are used which are valid at most one-half a period 
away. This enables the program to update the elements in steps of 2w, 
which results in a large time saving. The errors introduced by not 
completely updating the orbital elements are far less than the accuracy 
desired. 

(Hi) The determination of visibility at a site is not done by the 
obvious method of computing elevation and checking for a positive 
angle. The reason for this is that once the satellite's coordinates have 
been obtained, this method requires at least a square root, an arc tan- 
gent, and approximately thirteen multiplications. The method used 
instead requires only seven multiplications with the resultant saving 
in time. Instead of computing elevation, the program passes a plane 
through the site tangent to the earth. This plane, which can be con- 
sidered a ground plane at the site, has the equation 

aX. + /3F„ + yZ a - R„ = 
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where X„ , Y 8 , Z„ are the inertia! coordinates of the site, and R e is the 
distance from the center of the earth to the site. If, however, the coordi- 
nates of the satellite arc substituted in the equation, then a value other 
than zero is usually obtained. If this value is minus, then the satellite 
lies on the same side of plane as the center of the earth and is therefore 
below the horizon, but if the value is positive, then the satellite is on 
the opposite side of the plane which is above the horizon and is there- 
fore visible. The determination of visibility is thus reduced to evaluating 
a and /3, which are time varying due to the rotation of the earth, evalu- 
ating the four-term expression, and testing the sign of the result. The 
evaluation of a and /3 can be done using the previous values of a and /3 
with only four multiplications. It should again be noted that the price 
for this increase of speed is the loss of some accuracy. This method does 
not fully take into account the earth's oblateness. The result is that the 
plane is not exactly tangent to the earth, and a small amount of inaccur- 
acy is to be expected. (About 140 feet of error in placing the site.) This 
method is also limited in that it cannot predict rise and set at any 
angle other than 0°, but this information can always be obtained later 
from the COKE program. 

When the program has determined a period of visibility, it updates 
the orbital elements to either the center of the mutual visibility if any 
exists or otherwise simply the center of visibility, and punches out the 
elements plus other pertinent data on the pass. The program also prints 
out the pass number, rise and set times at the ground station, start and 
end of mutual visibility, and the maximum elevation seen at the site. 

This procedure is continued until a final time is reached. At this point 
the program finishes any pass it may be working on and then stops. 

9.2 COKE 

The COKE program uses the MUVIS results to generate an ephemeris 
that is exact but omits physical effects such as refraction and antenna 
distortion. Thus the tape can be generated ahead of time and just before 
the pass corrected for both meteorological conditions and the boresight 
corrections of the antenna to be used (there is only one antenna at 
each site at present). 

The program uses an Encke type method (see Ref. 0, p. 170) to solve 
for the satellite's position at four-second intervals. These are computed 
both forward and backward from the center of the pass or the center 
of mutual visibility, whichever the orbital elements have been updated 
to. Thus some rearrangement of data must be done to output the 
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ephemeris in time order. The heart of the program is the same inte- 
gration program used in the trend-removal portion of TELETRACK. 
Therefore only the peripheral programming needed to convert the results 
to pointing angles, to rearrange and output the results, and to control 
the direction and length of integration had to be written from scratch. 

9.3 AC EX PS 

The ephemeris tape required as input to the antenna digital control 
during a pass is generated by the expander program, ACEXP. In addi- 
tion to generating this tape, this program also produces the "mission 
printout," a listing of pertinent data regarding a pass over the ground 
station. The main input to this program is the data tape from COKE, 
containing time, azimuth, elevation and range at four-second intervals. 

One data point on the ephemeris tape contains the following infor- 
mation : 

(a) time 

(b) azimuth and elevation positions 

(c) azimuth and elevation first differences 

(d) azimuth and elevation second differences 

(e) azimuth and elevation predistortions, and 

(f) gain factor. 

Azimuth and elevation first differences for the *'th data point are com- 
puted according to 

where P, represents azimuth or elevation position. The second differ- 
ences are computed according to 

D 2 = P i+1 - 2Pi + P w . 

Azimuth and elevation predistortions, which are discussed in Ref. 7, 
are estimated to be functions of elevation only and to be of the form 

PD A = (a + bE)/cos E 

PD B = c + dE. 

Current values for the parameters are 
a = 0.015 degree 
b = -0.000786 
c = 0.057 degree 
d = -0.000712. 



1378 THE BELL SYSTEM TECHNICAL JOURNAL, JULY 1963 

The ground station transmitter gain factor, discussed more fully in 
Ref. 8, is computed as a function of range as follows 

7 -]|g«D +0.0791 (a- 1)] 

where a = log, S/S min for S mi „ < S < 10 S min 

= for S ^ S mia 

= 1 for 10 S m m < S. 

Smin is chosen according to the characteristics of the pass over the site. 
Elevations are corrected for refraction as follows. Index of refraction 
is computed according to 

n - 1 = (0.776 X 10~ 4 ? ; + 0.372 e/T)/T 

where T is temperature in degrees K, p is air pressure in millibars, and 
e is water vapor pressure in millibars (TCef. fl, pp. 13-15). The correction 

AE = (n - 1) cotE 

is added to the elevation, E, before putting it on the ephemeris tape. 

The mission printout is generated to aid the operating personnel 
during a satellite pass. Tabular data at one-minute intervals specify 
time, azimuth, elevation, range, one-second increments in azimuth, 
elevation and range, and Doppler shift. From a knowledge of the azi- 
muth rates the program predicts when (if at all) the horn antenna will 
lose autotrack due to excessive azimuth rates. The angular distances 
between the satellite and the sun are also computed, and if they come 
within 2° of each other an appropriate warning is included in the mis- 
sion printout. 

X. OPERATIONAL KKSULTS 

These programs have been a part of the Bell System satellite com- 
munications ground station operational system since the July 10, 1962, 
launch. Initial predictions were based on the launch and injection 
data, corrected by the few observations possible in the first six orbits. 
From the sixth orbit on, predictions were based entirely on track data 
acquired at the Andover site. By the seventeenth orbit (the second day), 
the orbital elements had been refined sufficiently so that the horn- 
reflector antenna autotrack could acquire the satellite using the predicted 
angles. From that point on the normal mode of acquisition was from the 
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predicted angles, and the use of the auxiliary antennas as acquisition aids 
was generally not required. This means that predictions have generally 
been within ±0.2°, at least at the horizon where acquisition is usually 

achieved. 

The launch of the Tclstar satellite was carefully planned to put the 
apogee in the northern hemisphere to maximize the periods of mutual 
visibility in the early phases of the experiment. During the first few 
weeks following the launch, the prediction accuracy was very good. 
Samples of the results of orbit determination and prediction during this 
period are shown in Table VII. The predicted angles, extending five 
days ahead, were generated from orbital elements computed from 
precision angular tracking data obtained during the preceding five days. 
The observed angles were obtained from the precision tracker. It should 
be noted that errors in azimuth should be multiplied by the cosine of 
the elevation in order to convert them to errors in sightline angle on a 
par with the errors in elevation. 

With apogee in the northern hemisphere, the tracking periods were 
long (over 30 minutes). As perigee precessed toward the northern 
hemisphere and the tracking periods became shorter, a gradual degrada- 
tion in prediction accuracy was noted. While the prediction accuracies 
were sufficient for the daily antenna pointing operations at Andover, 
they proved inadequate for providing pointing information for the 
optical experiment at Holmdel 10 and for determining satellite positions 
for the radiation effects study." These uses of the predictions require 
accuracies of 0.1° and both require that the satellite positions be related 
to geographical sites other than that at which the track data are acquired. 
This prompted a renewed study of the orbit determination method 
and the program implementation. This investigation revealed that this 
method is quite sensitive to observational bias, particularly when the 
track data are obtained from short passes rising to high elevations. This 
sensitivity can be reduced by using only tracking passes of 30 minutes 
or more in which the maximum elevations do not exceed 50°. However, 
that is a severe restriction to place on a single tracking site with a highly 
eccentric orbit such as that under consideration here. In addition, it 
was found that the approximate methods used to account for the per- 
turbations due to the earth's oblateness were inadequate except when 
the line of apsides is nearly parallel (as in July, 1902) or nearly perpen- 
dicular to the line of the nodes. Programs providing more complete 
perturbation calculations have been written and are presently under- 
going tests. 

From this study it was concluded that to achieve prediction accura- 
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Azimuth 


Elevation 


Date and 
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me 


(degrees) 


(degrees) 


Pass Number 


(hrs-m 


in, UT) 












Observed 


Predicted 


Observed 


Predicted 


7/23 


21 


35 


195.16 


195.20 


20.50 


20.48 


#124 


21 


41 


188.83 


188.86 


29.37 


29.38 




21 


47 


179.58 


179.60 


38.28 


38.29 




21 


53 


165.12 


165.14 


46.35 


46.37 




21 


59 


142.81 


142.80 


51.28 


51.29 




22 


05 


115.18 


115.21 


49.24 


49.22 




22 


11 


92.13 


92.15 


39.03 


39.04 




22 


17 


77.19 


77.19 


24.08 


24.08 


7/24 


21 


17 


187.80 


187.81 


25.73 


25.72 


#133 


21 


23 


179.40 


179.40 


34.40 


34.40 




21 


29 


166.95 


166.97 


42.47 


42.44 




21 


35 


148.30 


148.32 


48.26 


48.26 




21 


41 


123.75 


123.75 


48.91 


48.88 




21 


47 


99.92 


99.93 


42.22 


42.20 




21 


53 


82.60 


82.60 


29.77 


29.77 


7/25 


23 


43 


238.23 


238.24 


22.01 


22.03 


#143 


23 


49 


237.62 


237.61 


32.26 


32.29 




23 


55 


235.99 


236.00 


44.29 


44.29 


7/2G 


00 


01 


231.28 


231.31 


58.89 


58.94 


#143 


00 


07 


207.91 


207.83 


76.43 


76.44 




00 


13 


99.76 


99.82 


72.91 


72.91 




00 


19 


80.65 


80.65 


46.23 


46.22 


7/2G 


20 


31 


183.69 


183.68 


21.11 


21.12 


#151 


20 


37 


175.55 


175.59 


29.39 


29.42 




20 


43 


164.51 


164.54 


37.05 


37.04 




20 


49 


149.13 


149.17 


42.90 


42.91 




20 


55 


129.16 


129.21 


45.17 


45.18 




21 


01 


107.72 


107.77 


41.98 


42.02 




21 


07 


89.57 


89.59 


33.43 


33.44 




21 


13 


76.29 


76.29 


21.36 


21.35 


7/27 


20 


08 


181.81 


181.84 


18.66 


18.67 


#160 


20 


14 


173.80 


173.85 


26.84 


26.84 




20 


20 


163.25 


163.30 


34.25 


34.31 




20 


26 


149.03 


149.11 


40.15 


40.19 




20 


32 


130.89 


130.95 


42.99 


42.98 




20 


38 


110.85 


110.93 


41.09 


41.09 




20 


44 


92.90 


92.93 


34.18 


34.20 




20 


50 


78.98 


79.01 


23.62 


23.63 



cies of 0.1° or better, the angular observations must be taken from more 
than one geographical point or, if from a single tracking site, the angu- 
lar observations must be supplemented by an additional independent 
track measurement, such as slant range to the satellite. A program 
system avoiding the shortcomings of the present method is now under 
active development. This system uses a modified method of combining 
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passes and improved perturbation calculations, and has the ability of 
including slant range measurements and data from several tracker 
sites. 

The orbit determination method described meets the objective of 
minimizing the computer requirements by eliminating the mass storage 
requirements and time-consuming iterative procedures inherent in the 
classical differential corrections technique. As described, the method 
and programs are adequate for providing acquisition information for 
autotracking communications antennas if the tracking restrictions can 
be met. For a single tracking site, these restrictions imply a perigee of 
1000 nautical miles or more. If lower orbits must be handled or greater 
accuracies are required, the improvements mentioned above should be 
considered. 

APPENDIX A 

Derivation of Equations (7) and (8) 

Since /l _1 .f = (A' 1 + B~ x )x - B~ l x, (5) may be written in the form 
I = x - (A- 1 + B- 1 )- 1 B-\z - y). 



Now, 



Hence, 



(A -i + £-*)-' BT l = [B(A~ l + BT 1 )]- 
- [1 + BA" 1 ]' 1 

= [(A + B)A~T 1 
= A(A + B)-\ 

2 = x - A(A +-B)" 1 (x - y). 



Since, by (5), we may interchange x and y provided that .4 and B are 
also interchanged, we have 

z = y + B(A + B)~ l (x - y). 

Thus, if Wi and vh are any two six-by-six matrices whose sum is a unity 
matrix, 

z = Uh x + w$ - (u\A - w 2 B) (A + B)" 1 (x - y). 
Substituting A = PG~ l and B = QG~~\ and noting that 
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(A + B)- 1 = G(P + Q)- 1 



WiA - iv,B = (wiP - w 2 Q) G' 1 

we get (7). 

Noting that the right-hand member of (6) is a half of that of (5) 
if we replace x by A and y by B, (8) follows from (7). 
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